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INTRODUCTION 


During the past five years, numerous pioneering archival publications 
have appeared that have presented computer solutions of the mass-weighted, 
time-averaged Navier-Stokes equations (Farve, 1965) for transonic problems 
pertinent to the aircraft industry. These solutions have been pathfinders 
of developments that could evolve into a major new technological capability, 
namely the computational Navier-Stokes technology, for the aircraft industry. 
So far these simulations have demonstrated that computational techniques, 
and computer capabilities have advanced to the point where it is possible to 
solve forms of the Navier-Stokes equations for transonic research problems. 

At present there are two major shortcomings of the technology: limited com- 

puter speed and memory, and difficulties in turbulence modelling and in com- 
putation of complex three-dimensional geometries. These limitations and 
difficulties are the pacing items of the continuing developments, although 
the one item that will most likely turn out to be the most crucial to the 
progress of this technology is turbulence modelling. The objective of this 
presentation is to discuss the state of the art of this technology and suggest 
possible future areas of research. 

At present, the viscous transonic flow research is conducted by either 
a zonal viscous-inviscid interaction procedure or a global Navier-Stokes pro- 
cedure. There is no formal presentation of the state of the art dealing with 
viscous-inviscid interaction procedures at this Symposium. For this, one is 
referred to the proceedings of an AGARD Symposium on "Computation of Viscous- 
inviscid Iterations" (1980). These procedures have achieved some success 
but most either predict poorly or fail when faced with flow separation. The 
procedure of Le Balleur (1980) for small separated regions is promising. 

There does not appear to be a single one of these procedures which gives 
acceptable results under a wide range of conditions. Of course, these pro- 
cedures are being further developed, and in those cases where they can be 
trusted, they should be computationally cheaper to use than a global Navier- 
Stokes calculation. One expects that both the viscous-inviscid interaction 
procedures and the global Navier-Stokes approach will contribute to the 
understanding of various transonic flow phenomena and in providing insight 
for developing efficient numerical methods. 

We now discuss some of the flow conditions for which the Navier-Stokes 
equations appear to be required. On an airfoil there are four different 
types of interaction of a shock wave with a boundary layer: (a) shock- 

boundary-layer interaction with no separation, (b) shock- induced turbulent 
separation with immediate reattachment (we refer to this as a shock- induced 
separation bubble) , (c) shock- induced turbulent separation without reattach- 
ment, and (d) shock-induced separation bubble with trailing edge separation. 
The shock-induced separation is caused by a strong shock wave. A proper 
treatment of interaction of this shock with a boundary layer requires the 
Navier-Stokes equations, at least locally (Melnik, 1980). 

Shock waves that terminate in the vicinity of boundary layers are seldom 
steady, particularly on transonic wings and control surfaces. In some cases. 
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the shock-boundary-layer interactions are observed to oscillate periodically 
with relatively large amplitudes (Finke, 1975). These fluctuations can cause 
stalling, buffeting, flutter, and control-surface buzz. The first two phenom- 
ena arise at large angles of attack when the upper-surface separation of the 
boundary layer extends from the shock wave to the trailing edge and beyond. 

The last two phenomena are manifested when the separated boundary layer 
experiences lateral oscillations in the wake. A different type of transonic 
flow problem is recently reported by McCroskey, et al. (1981). They report 
transonic flow near the leading edge for free-stream Mach numbers as low as 
0.2 on an oscillating airfoil. This flow is characterized by a small super- 
sonic bubble with or without shock waves. At Mach numbers between 0.3 and 
0.5, the airfoil may experience shock induced leading edge stall (McCroskey 
et al. , 1981) . 

There are at least two motivations for understanding separated flows: 

(a) controlling and minimizing the effects of separation when it is an un- 
desirable feature, and (b) organizing separation so that it constitutes a 
natural way of improving aerodynamic performance. The latter occurs in 
three dimensions where strakes are used to create streanwise vortices that 
increase performance at cruise and climb conditions. It appears that air- 
craft designers are not so much worried about incipient or microscopic 
separation bubbles of small extent as they are about a boundary layer failing 
to reattach before the trailing edge. If that happens, it may cause, depend- 
ing on its severity, stall and buffet, pitchup motion, and possibly degrada- 
tion of lateral stability. 

When the boundary-layer assumptions are almost valid through a small 
separated region which is not caused by a shock wave, it is possible to 
determine, using the boundary-layer equations, the main effects of the sep- 
aration with an integral method (Le Balleur, 1980), and the quantitative 
structure of the separated region with a differential method. But when the 
separation region is not small, this approach fails, and the Navier-Stokes 
equations are required. 

In computational aerodynamics, both the physics and numerics are equally 
important. Physics is involved in selecting the appropriate governing equa- 
tions and formulating suitable initial and boundary conditions. Numerics, on 
the other hand, deals with generating a grid system, devising stable, accurate 
and efficient approximating schemes for solving the differential equations 
along with the initial and boundary conditions, and actually carrying out the 
solution procedure. All of the processes are important, and they all affect 
the accuracy of the solution. For the purposes being discussed here, the 
accuracy required of the solution is determined by the practical requirements 
of. the aircraft industry. If this solution fulfills these requirements, then 
it is accurate enough. The above processes dealing with physics and numerics 
for the Navier-Stokes equations constitute the Navier-Stokes technology. 

At present, computer simulations of transonic flow fields are usually 
validated by comparison with experiments which are in themselves simulations. 
This reliance on experiment results principally from the fact that the ef ects 
of turbulence must be modelled and the models are essentia y emp r ca . n 
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addition, this reliance results when a numerical solution, for all practical 
purposes, is not shown to be independent of the discretization errors. It is 
usually not possible to show the extent to which a large scale, numerical 
simulation is affected by discretization errors which is caused, at present, 
by lack of computer speed and memory. On the other hand, the validity of an 
experimental simulation is, more often than not, questionable. Generally, a 
quantitative assessment of effects of any known deficiencies in the data is 
lacking. Rarely are the initial and boundary conditions completely documented. 
There is usually a minimum rather than a comprehensive set of data. 

Keeping in mind the above general shortcomings of both numerical and 
experimental simulations, we discuss the state of the art of predictive 
Navier-Stokes technology dealing with the above processes and present some 
computed simulations of transonic flows. 


GOVERNING EQUATIONS 


Navier-Stokes Equations 

The continuum, compressible fluid mechanics is described by the class- 
ical Navier-Stokes equations, properly modified to take into account varia- 
tions in density and temperature, along with equations governing conservation 
of mass and energy and an equation of state, taken from equilibrium thermo- 
dynamics. This system is referred to here simply as the Navier-Stokes equa- 
tions. We shall assume that solutions of this system, subject to appropriate 
initial and boundary conditions, do exist and are unique. However, only local 
existence theorems in two- and three-dimensional problems have been estab- 
lished (Solonnikov and Kazhikhov, 1981); and the Cauchy problem for a perfect 
polytropic gas in three— dimensions is solvable "in the large" provided the 
initial data are close to constants (Matsumura and Nishida, 1980) . In short, 
the mathematical analysis of the above system is far from complete. 

In the Navier-Stokes equations, the assumptions concerning the stress 
tensor and the heat-flux vector exclude rarefaction shocks without specif- 
ically assuming the second law of thermodynamics. Therefore, the entropy 
condition (Lax, 1973) need not be satisfied by a numerical method for these 
equations. The effect of viscosity and heat conductivity develops a con- 
tinuous transition through a shock wave. In the transonic flow regime, 
these equations are valid through this wave which is, however, quite thin 
if its intensity is strong enough. For example, at a Mach number of 1.05 
and Reynolds number of 10^, the shock thickness in air is almost the same 
as the thickness of the linear sublayer of a turbulent boundary layer on a 
smooth flat plate. The latter thickness corresponds to about y'^~ 5, 
where y"*" is the Reynolds number based on the friction velocity and a 
length scale of turbulence. At lower Mach numbers the shock is even thicker. 
A shock wave with such a small thickness is not usually resolved in current 
transonic simulations. (Likewise, the contact discontinuity is not resolved.) 
Instead, it is considered to be a discontinuity, the location of which is 
part of the solution procedure. However, its thickness may not be small when 
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It begins to interact with a viscous boundary layer and can even lose its 
Identity as it penetrates into the viscous region. 


Reynolds-Averaged Navier-Stokes Equations 
with Mass-Weighted Variables 


In the study of turbulence by means of the Navier-Stokes equations, it 
IS usual to use some form of averaging. For example. Monin and Yaglom (1971) 

present a general space-time averaging procedure for functions f(x,t) given 
by the equation ' ^ B->-ven 


<f(x.t)> 


<X> 



f(x - c.t - x) g (? ,x)dc dx 


( 1 ) 


Here, the overbar and the underscore indicate an instantaneous value 
vector field, respectively. The non-negative weighting function, g, 
the normalizing condition 


and a 
satisfies 



g(l.x) d^d X = 1 


( 2 ) 


The choice of this weighting function determines the significance of the 
averaged quantities. For example, if g is a constant over some time 
interval T and zero outside of it, and the dependence on c is a Dirac 
delta-faction, <f(x,t)> is referred to a time-averaged quantity. In un- 
s eady flows, the interval T must be large compared to the periods char- 
acteristic of time scales that cannot be resolved computationally, but small 
compared to the period of resolvable flow motion. 


The system of equations determined by applying the above time-averaging 
procedure constrained with the Reynolds conditions (Monin and Yaglom, 1971) 

averaged Navier-Stokes equations. For compressible 
tiuids, these equations contain second-order moments, such as <p'u'> and 

<p'u’u’>, due to fluctuations iu the fluid density 
( an Driest, 1951). Here, the prime denotes fluctuating quantity. Therefore 
or these fluids, instead of time- averaged flow quantities, mass-weighted 
t me averaged quantities are preferable. For example, the mass-weighted 
velocity U£ equals to <pui>/<p>. This averaging procedure eliminates the 
above moments from the averaged Navier-Stokes equations but it does not 
remove density fluctuations from turbulence. This procedure appears to be 
irst used in the study of atmospheric turbulence by Hesselberg (1926) 

(Favre, 1969). A comprehensive discussion of this procedure for compressible 
turbulent flows is presented by Favre (1969) and by Cebeci and Smith (1974) . 
Henceforth, the equations resulting from this type of averaging are simply 
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called the Reynolds-averaged Navier-Stokes equations. These equations, with- 
out external forces, may be written in dimensional form as 


P + = 0 


(3) 

• 

(pui) + (pu.u^) . = -p . + 


(4) 

(ph) + (pu h) = p + u p + 

J J >J 

“i.j ■ 

(5) 

Here, p and p are, respectively, time-averaged mass density and pressure; 
uj^ and h are, respectively, mass-weighted mean velocity and enthalpy. The 
Cartesian-tensor summation convection is used. The overdot indicates a 
partial derivative with respect to time; and subscripts after commas denote 
partial differentiation. Further, the symbols and q j , respectively, 

represent the specific time-averaged total shear stress and heat flux as 
follows: 

‘’ij ■ 

-'ij 

(6) 

<pu'h' > 

, - ^ h 1 J 


(7) 

q . - n , 1 



where v is the kinematic viscosity. These include contributions of both 
the molecular and turbulent transport. The mean strain- rate tensor 
and the Reynolds stress tensor - pR-ij are given by 


S 


ij 





( 8 ) 



<pu^Uj> 

p 


(9) 


The above equations are identical to the equations used to determine 
laminar flows, except for the Reynolds stress tensor and turbulent heat 
flux vector, equations (7) and (9). In addition, these equations essentially 
exhibit a term by term correspondence with those for the incompressible 
fluids. This correspondence permits extension of the large body of expe- 
rience existing with modelling turbulence for constant-density flows to 
transonic flows, provided turbulence structure in both these flows is 
closely the same. Interpretation of Morkovin's hypothesis (Morkovin, 1964) 
suggests that this is the case for boundary layers and wakes at free-stream 
Mach numbers less than about 5 and of jets at Mach numbers less than about 
1.5 (Bradshaw, 1977). This hypothesis states that the effects of density 
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fluctuations on turbulence are small when the root~me an- square density fluc- 
tuation is small compared with the absolute density. Transport-equation tur- 
bulence models, which are discussed below in terms of mass-weighted, time- 
averaged variables, contain additional terms due to compressibility effects. 
These terms are negligible according to the above hypothesis in the transonic 
regime. 


Turbulence Modelling 

There are two approaches for turbulence modelling: the first-order 

approach in which the Reynolds stress tensor is modelled, and the second- 
order approach in which this tensor is determined from the Navier-Stokes 
equations. In the former approach, one forms the equations for the first- 
order quantities, such as mean velocities, and models the second-order quan- 
tities that appear in them. See equations (4) and (10). In the latter 
approach, equations are formed for the first- and second-order quantities 
(uj^ and Rii)> and the third-order terms are modelled. These equations may 
be simplified to yield algebraic stress models, which still require differ- 
ential equations, both for the turbulent kinetic energy and energy dissipa- 
tion (Rodi, 1980). 

In transonic, turbulent- flow simulations, the first-order approach is 
almost always used, and it forms the basis for the so-called zero-equation 
(algebraic), one-equation, and two-equation models. In practice, the actual 
form of these models and the manner of applying them generally differ in 
detail from investigator to investigator. General definitions and character- 
istics of these models are available from Cebeci and Smith (1974) , Reynolds 
(1976), Reynolds and Cebeci (1976), Rubesin (1977), Launder (1980), and 
Rodi (1981). (Simulation of transition is not considered in this state-of- 
the-art review . ) 

Some zero-equation models are based on the Prandtl mixing length hypoth- 
esis. But other first-order turbulence models are based on the "Newtonian" 
assumption, and they are, therefore, eddy-viscosity models. Boussinesq's 
eddy viscosity concept (1877) is based on an analogy with the gradient- 
diffusion mechanism of the kinetic theory of gases. Methods based on this 
concept are also known as eddy-dif fusivity or gradient-transport methods. 
Corrsin (1974) has presented limitations of gradient-transport models. In 
these methods, the eddy viscosity, v.j., is assumed to be a scalar and is 
defined by a Newtonian constitutive equation of the form 

\j = 2 - 3 (10) 

Here, v^ = is the turbulent kinetic energy. The v^ term may be 

absorbed in p. This relation restricts Rj^j and to the same prin- 

cipal axes, which is not true in general. It is possible to modify this 
relation in order to remove this restriction (Saffman, 1974). Algebraic 
models relate vj directly to Reynolds-averaged field quantities. Both 
one- and two-equation models contain a partial differential equation 
for turbulent kinetic energy, which defines a turbulence velocity scale, v. 
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One-equation models use a prescribed, empirical length-scale distribution, 
and two— equation models use an additional partial differential equation to 
define a turbulence length scale, Jl. A combination of the turbulence 
velocity and length scale determines the value of the eddy viscosity 


= evil 


( 11 ) 


where c is a constant. Methods using partial differential equations for 
turbulent quantities are also called transport-equation methods. 

In the eddy-conductivity concept, the transport of heat due to the 
time-averaged product of fluctuating enthalpy and fluctuating velocity is 
modelled. It is assumed that the turbulent heat flux follows a law similar 
to Fourier's law. Further, it is generally assumed that dynamic eddy vis- 
cosity, M™, and turbulent thermal conductivity have the same functional 
relationship with temperature. Although the turbulent Prandtl number varies 
across the boundary layer, it is commonly considered to be a constant, and 
it is usually taken to be 0.9 for air. Apparently, more complex modelling 
of the turbulent heat flux than this has yet to be attempted in transonic 
simulations. 

There are many zero— equation models. As an example, a model used by 
Baldwin and Lomax (1978) for attached, separated, and wake flows is briefly 
outlined below. This model is patterned after that of Cebeci (1971) . The 
turbulent boundary layer is regarded as a composite layer consisting of 
inner and outer regions. In each region, the distributions of v and 8, 
are prescribed by two different empirical expressions. For example, in the 
log-law region, £ is proportional to y, the distance normal to the wall, 
and in the outer layer, £ is proportional to the boundary-layer thickness. 
The proportionality of £ to y is extended into the viscous sublayer with 
a damping function suggested by Van Driest (1956). In the outer region, 
the vorticity is used to define the boundary-layer thickness. 


In the inner layer, 0 5 y 5 y^,* the expressions for v and £ are 


(v) 


inner 


£ 1^1 


( 12 ) 


and 


inner = *^1^ly26v^) ] 

with c = 1.0 in equation (11). Here, = 0.4, Q is the vorticity, and 
subscript w indicates wall values. 

In the outer region, y > y^. the expressions for v and £ used by 
the Baldwin-Lomax (B-L) model are: 
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(v) 


outer 


-•inax 
or 

0.25 U^,,/L 

dir tnax^ 


the 

smaller 


(14) 


and 


(J?-) = y C 

^ outer ■’^max BL 2 


(15) 


with c the Clauser constant equal to 0.0168 in equation (11). The quan- 
tity is the difference between maximum and minimum absolute velocity, 

the value of CgL is 1.6, and the Klebanoff intermittency factor, a 2 , is 
given by 


a 


2 



(16) 


The quantities y and L are determined from 
max max 


L(y) = y|^2|[l - exp (-y/ja.„| /26v ) ] (17) 

' 12 'w w 

The above exponential term is negligible in the outer part of the boundary 
layer. In wakes, it is set to zero. The quantity L^ax is the maximum 
value of L(y) that occurs in this equation, and y^ax is the value of y 
at which it occurs. 

The region of validity of the inner and outer scales is determined by 
y^* It is the smallest value of y at which values of inner and outer eddy 
viscosity are the same. The value of in the inner region and of c in 

the outer region are assumed to be universal constants for Rg > 5000, where 
R 0 is based on the momentum thickness. At lower Reynolds numbers, they are 
functions of Reynolds number. 

As an example of a two-equation model, the Wilcox-Rubesin (W-R) model 
(Wilcox and Rubesin, 1980) is presented below. This model is an extension 
of the model developed by Wilcox and Traci (1976) , which evolved from the 
model formulated by Saffman and Wilcox (1974) and that by Saffman (1970). 

In the earlier models, the term determining the rate of production of kinetic 
energy was inconsistent with that in a stress-equation formulation. The 
present model removes this inconsistency. In this model, the turbulent 
kinetic energy and the specific energy dissipation are given by 

(pv 2 ) + (pujv 2 )^^ = 2 pa^jU^^j - B^po)v 2 + [(y + 62 '^T^'^O ^ , j 
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(19) 


• 2 . 


2 


+ [ ( y + g^Vrp)^ ^ ] 


>3 >3 


where the length scale is defined by 


i = 


V 

OJ 


( 20 ) 


The eddy viscosity is computed from equation (11) and the constitutive equa- 
tion (10) is used to provide R^j • Wilcox and Rubesin (1980) recommend 
following values of the constant in the above model: 

6^ = 0.09, ^2=33 = 0.5, 8^ = 0.15, 3^ = 3^ = ^ 

C = [1 _ (1 - 3^) exp (- Re /2) ]/2 
D I 

6^ = 3^ [1 - (1 - 6^) exp (- Re^/4)]/c 


The turbulence Reynolds number is calculated as 


Re 


T 


y_e 

V 


The Reynolds-averaged Navier-S tokes equations along with turbulence-model 
equations constitute the governing equations of the Navier-S tokes technology. 


Conservation-Law Forms 

The Reynolds-averaged Navier-S tokes equations, as presented in equations 
(3) to (5), are not in a form generally suitable for simulations of flow fields 
around aerodynamic shapes. For such shapes, surface-oriented coordinates are 
preferred. Furthermore, the choice of dependent variables made for these 
equations is not the only choice available. 

For unsteady flows, Moretti (1979) recommends using the velocity com- 
ponents, pressure (actually In p) and entropy for the dependent variables. 

This is motivated by the fact that, in inviscid flows, there are two types 
of surfaces across which flow quantities can be discontinuous; character- 
istic surfaces and stream surfaces. Across the characteristic surfaces, 
pressure and velocities are discontinuous, but not entropy. In contrast, 
across stream surfaces, entropy is discontinuous, but not pressure. Any 
other thermodynamic parameters, such as energy or density, are discontinuous 
across both the surfaces. The above recommendation does not, however, lead 
to the conservative-law form (Lax, 1957 and 1973, and Richtmyer and Morton, 
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1967) of the governing eauations Tf 

variables may be appropriate. Anothfr cbni crucial, then the above 

the contravariant components of the velocitv^vec t use density, energy, and 
divergence form of the equations It is rn« leads to a non- 

the divergence or conservative-law form (v- these equations in 

Stone, 1980). A third choicTis to use Eisen^n and 

and conservative variables. It is this components of velocity 

conservative-law form that is presented below."" 

The conservative-law form of the Navier- 9 ^^l, 
variables facilitates capturing of discontinui^ ^ equations in conservative 
conservation of fluxes. The importance of r) naintenance of global 

cations and acceptable error boLds L Ld^car ^PP^i- 

of Navier-Stokes equations, shock waves and earlier discussion 

physical discontinuities during flow limnW regions are treated as 

form of these differential equations avoids above conservative 

discontinuities. Further therp ■{<= fictitious sources along these 

the absence of dif feren^ilbilit; acror^hrd"^”" equations i" 

these theoretical results facil^rarf discontinuities. In principle 

these theoretical advantages are maintained^"^ discontinuities. Whether’ 
upon the numerical scheme along with the grid^''^^ '''' ^ depends 
in the next two main sections of thL papfr LiL " "!! discussed 

taming global conservation of fluxeJ Hp^^^a Likewise, the issue of main- 
Analytical integration of a con^Ltive 1l r numerical scheme, 

pendent variable yields a difference be^i^ k’"" respect to an inde- 

Construction of a differencing scheL thiror^"'' 

erty is relatively easy if thf cons^a^ive ^ " conservation prop- 

Sometimes it is possible to formulate a differLci^'" is used to begin with, 
fluxes, starting with the nonconservative f^r^ Jo® ^ conserves 

It IS the discrete form that governs the rr, However, strictly speaking, 

differential form. Global consIr^^tioXrf 

assure that the discontinuities are ranr does not automatically 

are analogous to truncation errors TTont ’ ^on servation errors 

bounded and do not affect acceptable a«iracv^",- 
not the governing equations are in the conser^ative!!™".^^"^ 

to o. a.o ap.., 

Stokes equations.’ If ’a"""'’’' — 

law form, then similar numerical trer^m^n^ formulated in the conservation- 

equations of Navier-Stokes techLlogJ “ govemlnr 

formed'’from''thfSttef^^“oo?«^JS to '""“tvative-law form are trans- 
the resulting equations are not automatical / coordinates, 

although they can be made to be so (Viviand^^ conservative-law form, 

to facilitate global conservatL orfluxL’ 

form is not necessary for obtaining thp v ^^oretically , however, this 

for avoiding fictitious sources along + v.^) -> 0 and 

coefficients multiplying the transformed S:"ima"t:mm’mnSmhefr‘'f 





derivatives with respect to the new independent variables are continuous. 
This can be readily demonstrated following Lax (1954) , but the derivation is 
not given here. Further, it can be shown that the shock speed in the 
Cartesian coordinates and the curvilinear coordinates differ by a factor 
containing the metric coefficients. 


Reynolds-Averaged Navier-S tokes Equations 
in Curvilinear Coordinates 

Below, the Reynolds-averaged Navier-S tokes equations in conservative 
Cartesian variables are presented in arbitrary curvilinear coordinates 
In nondimens ional form, these equations can be written 


iQ 

8 T 


d 8C. 

T — 


Re 


d 3 V. 

T. — 

ki 


( 21 ) 


where d is the number of dimensions, and Q, C, and V are vectors 


Q = [p, pu^, ... , pu^, e] 


C. = QW. + p «>. 
1 1 1 


d 3 5. 

V. = ® L R- 

1 j 3x. 

J=1 J 


35 . d 3 5^ 

*'1 ■ IT * “j TT 

J=1 J 

^ '^jl’ ’ "^jd’ 

J=1 J 




O , . 


= (y + Z ( - T 


U r 

’ id" " 1 11 

J=1 


4- ( i V ■; 

i^iV ^ JTi ‘‘j 

3 5, 3u. 3 5, 3u.\ 

3 X . 3 5. 3 X . 3 5, / 

J k 1 k' 


( 22 ) 
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^1 = ~ M.) f '"I 


" Jtl ax 


^ r ^ 

*= tl 2 


P = (Y - l)p« 


a ^-r d 9 C . 3 X. 

— - - E — — L 

9x. 9 t 
J = 1 j 


9 t 


Q> = 


9 (x. , ... ^ X ) 


3 . 


^d> 


!!i 

3x . 


® «I+ 2 ' 


In the above expression, subscriots -i+i j /• . 

in a cyclic order, (1 2 3) (2 2 T» t- ’ tu^ J+2) vary 

(3A + 2U) of local ther^na^’ ^ hypothesis, 

per-unit ;»^Li and >■“= '•'=»" and total energy 

and ej. respectively energy-per-unlt sass are represented by e " 

lence stress ter„s. The .on,:„tu: :ru^ " f 

surface (? 2 ~direc tion) is retained ^ Tf ^ direction away from the 

ctss“c:^ b^fd’ app”-rt\o“r4S°ir\tio^^ja^ 

thln-shear-layer approxWlon^j^Uf J Yt^Z'the^basC th^^th''* 

gS\;LJua:ra;;:j:y:’":“hrke“; 

valid only for •'snail" separation bubbled and for 'vS^^sJock^bo ' 7 '''"^' 

approxlnatlon applied to e,uatlonl 2 ir?:ads“tr?S'f:L" 


12a •f i!i a ^ ’«2 

iti Hj Re 3x 3{, 


(23) 


with 


12 


T 

CT. . 

IJ 


Cw + 


_ 2 . 
3 h '^2 




Y 

Pr 




"'t/ 9x^ 3^2 


( 24 ) 


Instead of equation (23), some investigators (e.g., Steger, 1978, and 
Pullium and Steger, 1980) use, for convenience, the following equation 





(25) 


Boundary Conditions 

Boundary conditions for the above governing equations are determined 
by mathematics and physics. The mathematical character of these equations 
dictates the number and type of these conditions that determine the well- 
posedness of these equations. Further, this mathematical character is 
determined by the theory of characteristics. Theoretical analyses of two 
kinds are available: one based on the classical energy method (e.g., Elvius 
and Sundstrbm, 1973) which follows the earlier work of Serrin (1959), and the 
other on the normal mode concept (Kreiss, 1970). Most of the work is done 
for the compressible Eulerian and shallow-water equations, A few recent 
studies deal with the compressible Navier-Stokes equations (e.g., Oliger and 
Sundstrbm, 1978, and Gustafsson and Sundstrbm, 1978). These studies consider 
both number and a possible set of admissible forms of the boundary conditions. 
At present, such studies serve as a guide rather than as a useful tool in 
practical transonic simulations. A theoretical study of the well-posedness 
of the governing equations of the Navier-Stokes technology has yet to be 
done. The boundary conditions discussed below are based on both the math- 
ematical character of the equations and physical considerations. They are 
not based on the analytical procedures mentioned above. 


The mathematical character of the system represented by the linearized 
form of equation (21) is incompletely parabolic (Belov and Yanenko, 1971) or 
parabolic-hyperbolic. Without the time derivative, it is elliptic-hyperbolic 
^e system given by equation (23) or (25) is incompletely hyperbolic or 
hyperbolic parabolic. This system is parabolic only in (C? ~ t) plane. The 
global character of these systems remains the same even if the local character 
may be, for xnstance, purely hyperbolic. Therefore, the boundary conditions 
are determined by the global character of these systems. 


First we discuss the boundary conditions for the system represented by 
equation (21). Consider each equation of this system separately from the 
others as an equation determining Qj^; the other Q's in this equation are 
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assumed to be known quantities. Here Qi is a component of vector Q. The 

^rLtiorflf T condition in each coordinate 

direction for The second derivative of in Qo-equation requires 

conditions in each coordinate direction. Likewise, two boundary 
conditions are required for the remaining Q's. This means that if 3-3? is 

spLifyinr^O° region , then everywhere on conditions 

^ ® required; and on a part of 3.'^, a condition 

p c ying IS needed. These considerations determine the number of 

boundary conditions on . The type of the boundary condition for a Q. 
in any direction is determined by the highest derivative of this 0. in 
that direction. The boundary condition should be one order lower than the 

constraint yields boundary conditions which are 
either Dirichlet, Neumann, or mixed type. 

heuristic considerations help formulate boundary conditions 
based on physics. A set of these conditions for equation (21) are presenLd 

aries’ 3 « °'*a 1™ hmds of boundaries arise: rigid-wall bound- 

fCld’alo^e’ S "P™ ■ The rigid wall constrains the flow 

field along . This physical constraint is relatively easy to formulate 

and convert into computational boundary conditions. Open bounLries do not 

provide a material constraint, and hence appropriate conditions are not 
obvious. 


boundary provides velocity and temperature conditions on 

iS's’thai lS-2rr’' ^ ordinary conditions (Knudsen numbers 

less than 10 ) is accurately described by the no-slip and no-temperature 

jump conditions. These are the only two physical conditions avaiLb^e (In 

namelT ’ n only one physical boundary condition, 

namely, no flow normal to the rigid walls. Further, for an inviscid flow 

past an airfoil a Kutta condition must be imposed at the trailing edge of 

the caee of Impermeable walls, the no-flip Ln- 
Ffffhfr fr vanishing contravarlant velocity components,-*; - 0. 

The mass conservation equation governs the material derivative of 0-, 
Consequently, on 3.^, changes if its previous history is know^! otheSilse 
a condition on must be specified. This means that if fluid is on ' 

fluid^ent ' * determined by the mass conservation equation. But if 

fluid enters Jf by crossing 3.5P, must be specified. Therefore 0, 

fSfftlfa''’'Sb fb°" from its mfffhai 

r vative. When this recourse leads to numerical difficulties a new 

fqua^onf trioj^the appropriately combining the 'momentum 

. r derivative of pressure. After expressing 

pressure in terms of Q*s (equation of state), we have 
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The left-hand side of the above expression simplifies for orthogonal curvi- 
linear coordinates. As equation (26) is derived from the momentum equations, 
and as it replaces the mass conservation equation, it is not a boundary con- 
dition on Qi. This equation is subjected to the no-slip condition, when 
It is used. The above viscous terms vanish when the first-order thin-shear- 
layer approximation is valid; and they can be neglected only when the second- 
order thin-shear-layer approximation is valid. 


These conditions are also valid for internal flow problems. However, 
when simulations of the external flow problems include wind-tunnel wall 
effects, one alternative is to use the no-slip condition. Another alternative 
is not to compute the wall boundary layers. In this case, obviously walls 
cannot be considered as open boundaries if they interfere with the flow field 
around an aerodynamic body from that observed in free flight. This is the 
situation of present transonic wind tunnels. An ideal situation is to 
measure all required flow quantities just outside the wind-tunnel wall 
boundary layers and use these values as boundary conditions. Probably the 
next best avenue is to measure only pressure, again perhaps just outside the 
wall boundary layers, and then consider the boundary formed by pressure 
measurement locations as an open boundary. Another approach is to contour 
the wind-tunnel walls, such that they coincide with streamlines in free-flight 
conditions. The slip boundary condition is enforced along these contoured 
walls. This is restrictive, because in unsteady flows these free-flight 
streamlines, at a short distance from the body, can be time dependent. In- 
stead of these alternatives, the adaptive wind tunnels (see for Instance, 
Sears, 1981) could allow the use of the free flight boundary conditions.* 

The inflow, outflow, and tangent flow open boundaries require different 
treatments. The above discussion dealing with material derivative of 
shows that on inflow boundary, must be specified. On outflow, Q, is 

determined from the mass conservation equation; and on tangentflow boundary 
equation (26) is used. 


For external flow problems, boundary conditions are available at infin- 
ity, but not at finite distances. If the inflow boundary is at, say, about 
ten times the characteristic length of an aerodynamic body, then the influence 
of the body at that distance should be negligible, and therefore, it is 
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resolved frequencies beyond the error bounds for these frequencies. One 
point of view is that these terms do not change the mathematical character 
of the original equation because they are a part of the truncation errors of 
a numerical method. (See page 331 of Richtmyer and Morton, 1967.) These 
terms disappear when 0. An alternative point of view is to consider 

these terms as part of the original differential equations, since A? is 
never equal to zero. In this case, these terms do not change the character 
of the original parabolic equations. Also, they do not change the global 
character of the original hyperbolic equations, provided they do not intro- 
duce any boundary layers at the boundaries. This is achieved by not adding 
these terms either on the boundaries or next to the boundaries in ^2“ 
direction. This avoids additional boundary conditions for both parabolic and 
hyperbolic equations. These terms may form interior ’’boundary layers" such as 
captured (smeared) shocks. In this case, the "additional boundary conditions" 
for these terms are automatically provided by the appropriate neighboring, 
interior flow quantities. 

Some numerical methods require extra boundary conditions (e.g., Mehta, 
1977, and Yee, 1981), These condit ions are called numerical boundary condi- 
tions. 


COMPUTATIONAL GRIDS 


A computational grid system is a necessary part of any numerical solu- 
tion based on a finite difference, a finite volume, or a finite element method. 
The selection of a grid system is based primarily on the requirement for 
accuracy in the final solution. Secondary considerations are the effect on 
computational efficiency of the solution algorithm, and finally the ease of 
grid generation using available computer architecture. These concepts are 
discussed below. 


Accuracy Requirements 

Accuracy requirements are determined by the application of the numerical 
solutions of governing equations along with initial and boundary conditions. 

If the solutions serve the purpose for which they were intended, then the 
accuracy requirements are satisfied for that particular application. These 
requirements vary with purposes of applications and frequently tend to be 
subjective. Unlike the accuracy requirements, the discretization (truncation) 
errors are independent of both purposes of applications and subjectiveness. 
Therefore, in the discussion that follows, the accuracy constraints are not 
quantified, and the emphasis is placed on the discretization errors. 

Simulations of flow regions, throughout which the scales of motion are 
essentially the same in all directions, are probably best carried out by 
equi-spaced Cartesian meshes. In this case, the evaluation of mesh errors on 
the solution is completely determined by the size of the single-space inter- 
val. On the other hand, a flow field with a surface along which there is a 
turbulent boundary layer, is generally computed using a highly "stretched" 
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available, because computational resources do not permit it, A necessary but 
far from sufficient condition is to make sure that a given grid, along with 
a numerical scheme, maintains the free stream if the free-stream boundary 
conditions are applied. 


Methods For Generating Grid Systems 

At the present time, the property of providing a body-fitted grid system 
for external aerodynamic problems is automatically satisfied in two dimensions 
by all currently popular curvilinear grid-generating schemes. But appropriate 
interior grid systems for each class of topological geometries require trial- 
and-error manipulations of different variables in these schemes. (See, for 
instance, Sorenson, 1980.) Based on methodology, there are two types of grid 
generating schemes, algebraic and differential. An algebraic grid-generation 
scheme is a direct approach. It may be further classified into a conformal- 
transformation procedure and a nonconf ormal-transf omiation procedure. A 
description of conformal transformations for computational aerodynamics is 
given, for example, by Sells (1968), Ives (1976), and Morotti (1980). (A con- 
formal transformation may be defined either by an analytical function or by 
two Laplace equations resulting from the fact that the real and imaginary 
parts of an analytical function are harmonic. The procedure based on the 
latter definition does not require a separate discussion.) Some of the non- 
conf ormal procedures are the parametric multisurface transformations (Eiseman, 
1978 and 1979), transfinite interpolations (Eriksson, 1980), and the iso- 
parametric mappings (Forcey et al. , 1980). (Note, Eiseman has not used the 
adjective, "parametric.”) On the other hand, a differential grid-generation 
scheme is an indirect approach. This again may be further categorized as 
that based on a hyperbolic differential system and on an elliptic differential 
system. A hyperbolic procedure was first presented by Barfield (1970); and 
then it was extended and analyzed by Starius (1977). Recently, Steger and 
Chaussee (1981) have modified Starius* procedure. Thompson, Thames, and 
Mastin (1974) exposed the elliptic procedure to the computational aerodynamic 
community by extending, in particular, the work of Barfield (1970), Godunov 
and Prokopov (1972), and that of Amsden and Hirt (1973). 

When a boundary of a flow field can be mapped with an analytical function, 
when the resulting distribution of boundary grid points is nearly satisfactory, 
and when the interior grid distribution is less of a concern, conformal trans- 
formations are the best. They give rise to simple geometrical mapping quan- 
tities, and it is easy to assemble a grid system with them. Furthermore, 
they provide exact values of geometrical quantities. These transformations, 
however, cannot be extended to three dimensions, but they can be used in two- 
dimensional cross sections of a three-dimensional flow field. 

The hyperbolic transformation procedures give, in two dimensions, orthog- 
onal curvilinear grid systems. With these procedures, it is not automatically 
possible to control either the location of the outer boundary or the distribu- 
tion of points on it. Therefore, they cannot be used directly for internal 
flow problems or for patching different grid systems. Further, their applica- 
tion and usefulness in three dimensions remain to be demonstrated. On the 


19 



thrpp elliptic transformation procedures have been extended to 

three dimensions (e.g. , Mastin and Thompson, 1978, and Lee et al 198n^ 

tlonal time than her proLLres ^^ar“ r ”” 

all^i Wa? '>'? of grid system at the boundaries. But they do 001 *^ 

the :uwL*““„rrL”' system beLuse of 


tion p^o^L^e: auoS“:c:rg:Jr::;t^o!.^%T:^:‘:e™“at“r:^:f t‘h”“^°™- 

Of generating vLubles "ban ^rr^thers!’' =»-‘»-tlo„ 

^ Sometimes it -is possible to choose the type of grid pattern For 

oriH TK-ir. • - j grid. One may also use the 'H' 

?erticariLrir?^r^« ieooetrical singularity, if the two halves of a 

n the H, one below and one above the horizontal lin#=* moo<- 
at an angle other than 180° as in ' ’ TVn- o * line, meet 

“rifif^s 'j; of^sr^eirsLstaiiJy"^^:"'; 

tL ^egL-b JL^d^Jhf t-:?inret°r :If‘1^rea-: fL-hi ^1"-gr“: 

rtr=:n°gtg=f^:: r^'jo ;Ln ' 

'o' three grids as they are usually programmed the 

ar"a980rh'^ r’'"d LeranrSert^ufsoran^Le 

eL^^biocr^Tirtv 

. th this approach, there are two major shortcomings It Intro 
duces geometrical singularities In the transformed domain wh«e therfwere 
none to begin with and the grid control In the physleal dc^alris poor! 


^In contrast, current inviscld transonic 
with ’o' grids that use about 150 grid 


simulations are generally conducted 
points on the airfoil surface. 


20 



particularly, across block boundaries and along the trailing edge of the wing. 
The above investigators have also considered a single-block around a wing- 
body configuration. In this case, geometrical singularities become regular 
in the transformed domain. Eriksson (1980) has used an algebraic scheme for 
the same configuration. The resulting 'C grid pattern around both the 
leading edge of the wing and the wing tip appears to be acceptable. Moretti 
(1980) has shown how to assemble a grid system in cross-sectional planes of 
a fuselage-and-arrow-wing configuration using conformal transformations. 
Complex three-dimensional geometries are first rendered quasi two-dimensional, 
then two-dimensional grid-generating techniques are applied. 


Methods for Improving Flow Simulation Accuracy 

One requirement of accurate solutions is that they be, for all practical 
purpose, independent of the grid system. So far, this has not been system- 
atically demonstrated for turbulent, transonic simulations. The generally 
accepted practice of indicating the order of truncation error of a numerical 
method does not quantify the discretization errors. Although quantification 
of these errors is difficult, it is possible to determine their effects 
through grid— ref inement studies. On the other hand, minimization of these 
errors may be achieved by a proper choice of both the numerical method and the 

system. Usually, there is more freedom in choosing the grid system than 
in choosing the accuracy of the numerical method. Further, the choice of the 
grid system is determined by a priori knowledge about the solution. Most of 
this knowledge is available in terras of generalities rather than specifics. 

For example, surface boundary layers are always resolved with the help of 
some stretching function near the known surface. But without the specific 
information, such as the magnitude and location of gradients in the flow 
field, the grid system employed can often be wasteful and not satisfactorily 
concentrated on those regions where a better resolution is desirable. 

For a better utilization of grid-point resources, there is a growing 
interest in solution-adaptive grid systems. In a moving finite-element 
method, which allows both nodal amplitudes and nodal positions to move con- 
tinuously with time, nodes generally move automatically to those regions 
where they are most needed (Gelinas et al. , 1981) . In finite— difference 
methods, there are currently two basic strategies. The first strategy 
involves tracking a fluid property, such as the density gradient, and insert- 
ing or regridding so that finely spaced grid— points are in the immediate 

vicinity of that selected property (for instance Dwyer, 1980, and Kovenya and 

Yanenko, 1980). The second strategy is to minimize the leading term or 

terms of the modified equations that determine the order of the truncation 
error of a numerical method (e.g,, Pierson and Kutler, 1980, and Rai and 
Anderson, 1980). 

So far, the adaptive grid techniques have been primarily applied in one- 
and two-dimensional Burgers’ equation, and for a two-dimensional heat equa- 
tion. Extension of these techniques to the Navier-Stokes equations for 
turbulent, transonic simulations is a difficult undertaking. Questions, such 
as what flow variables to monitor, which truncation errors to minimize. 
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and zonal methods and finite element methods do not necessarily produce well 
ordered data bases. The real importance of well-ordered data bases occurs in 
studies involving three-dimensional spaces, and we have very little experience 

in this area. 


NUMERICAL METHODS 


Two Crucial Nonlinear Convective Phenomena 

In order to clarify the discussion presented below, it is useful to 
develop a concept that can be used to relate physical and numerical phenomena. 
We search for some form of scale in both time and space that is common to 
both phenomena, and find an excellent candidate in the 

harmonic analysis made of the physical variables with reference ^ther tim ^ 
or space. The physical side of this concept can range from the very natural 
(in Lperimental studies of isotropic turbulence) to the rather contrived 
(in the harmonic analysis of a discontinuity). On the numerical side, these 
Lequencies form part of the exact solution to certain model linear problems 
with periodic boundary conditions, but are only loosely related to the eigen- 
system of most difference equations actually being solved. Nevertheless, the 
association of frequency with scale is a very convenient concept when discuss- 
ing some of the broader aspects of the numerical simulation of fluid flow. 

The Euler equations model an unsteady flow that can contain a discon- 
tinuous solution referred to as a shock wave, or simply as a shock. For the 
Navier-Stokes equations, shock waves are not, strictly speaking, discontinuous, 
their thickness being of the same order as the thickness of the linear sub 
Lyer In a turbulent boundary layer (see the section. Navier-Stokes Equations), 
ile spectral analysis of a variable having a discontinuity, or an abrupt jump 
that is "nearly" discontinuous, is shown in figure 1. Notice that all, 

"nearly" all, of the high-frequency terms have finite amplitude. In the 
themrof the previous pLagraph, all or nearly all scales are P-sent. This 
has an important influence on the construction of numerical methods used to 
compute flows with embedded shocks. 

In this paper, we are interested in flows that have significant regions 
of turbulence and separation. Laminar flows and flows with attached tur u 
lent boundary layers can be computed using the methods we are discussing 
they usually can also be calculated by simpler and less expensive methods. 
Although the vorticity that is essential for the production of turbulence is 
generaLd by the viscous properties of the fluid-surface interface and curved 
shock wave, turbulence itself is generated away from the surface and caused 
by the nonlinear interactions of the convection terms in the Euler equations 
the same terms responsible for the generation of shock waves, ^or the point 
relative to this discussion, the most illuminating aspect of 

lies in the spectral representation of its inertial range sho^ in figure 2. 
This gives the amplitude of the kinetic energy associated with each harmonic 
in a spectral analysis of a typical high Reynolds number turbulent 
TennekL and Lumley (1972). Notice that the scales of both axes in the figu 
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are logarithmic, that almost all of the energy is carried in the low wave 
numbers, and that molecular dissipation is limited to the relatively high 
wave numbers where the energy content is low. The flow represented by the 
high energy, low-frequency region is referred to as large scale or large eddy 

represented by the opposite end as low scale or small 
eddy motion. The results shown in figure 2, referred to as the energy cas- 
cade, greatly influence models used to approximate the effects of turbulence. 

Numerical Techniques for Computing Shocks 

There are two common approaches used when devising numerical methods for 
calculating flow fields with shock waves. They are referred to as shock 
fitting and shock capturing. Shock-fitting methods employ some kind of test 
for detecting the shock location, and then treat the shock as a local dis- 
continuity across which the Rankine-Hugoniot relations must be satisfied. 

oc -fitting methods are probably to be preferred where they can be generated 
by reliable and efficient codes. They eliminate the need for conservation- 
law forms of the governing equations (which has certain simplifying attrac- 
tions), and they produce sharp discontinuities at the jump location. They are 
quite popular for computing many flows that can be modelled by the inviscid 
Euler equations, especially where the flow field is supersonic, see for 
instance, Kutler (1974) and De Neef and Moretti (1980). However, the flows 
of interest in this report can have strong shock boundary-layer interaction 
and the effect of viscosity must be included in this region. Further we are 
interested in the flows that contain three-dimensional and oblique shocks. 
Shock fitting under these conditions can become extremely difficult, and our 
remaining attention is limited to shock-capturing methods. 

The point of a shock-capturing technique is that the shock forms and 
moves about in a mesh, while some kind of. analytic connection is maintained 
between the flows on the two sides of the wave front. This does not mean 
that the shock-capturing methods cannot have built-in logical tests that try 
to isolate the shock location. Very often they do, and very often they make 
use of the test results to make local adjustments to the differencing scheme 
to improve its capturing capability. Still, by definition, a shock-capturing 

numerical method connects the dependent variables on the two sides of the 
wave. 

An immediate consequence of shock capturing relates to the spectral 
structure of a discontinous function shown in figure 1. Since the capturing 
technique is based on some kind of numerical continuity across the shock the 
harmonic analysis can be used to represent the result. It is well known ’that 
a inite grid can only support a finite number of frequencies in a discrete 
ourier series. For example, an equispaced grid of M points can accurately 
accommodate k - M/2 harmonics of the form eikx. Frequencies higher than 
k reappear as lower frequencies, a property referred to as aliasing. In an 
unsteady flow with a moving shock, these higher frequencies are constantly 
eing generated by the nonlinear convective interaction. For example the 
product of the waves e^^x giix brought about by terms such as u 9„’v 
produces two harmonics, one having a lower frequency proportional to k’- i. 
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and the other having a higher frequency proportional to k + Ji. This behav- 
ior can be verified in numerical simulations by observing how a simulated 
shock constantly tries to steepen. A linear discontinuity shows no such 
tendency. The situation just described can be summarized as follows: 

(1) Any discrete grid system can accurately support only a limited 
number of low frequencies. If higher frequencies are placed on it, they 
appear as amplitudes of low order terms. 

(2) Convective nonlinear interactions are constantly cascading low 
frequencies to higher ones. 

The numerical difficulty brought about by this situation in the case of 
shocks is illustrated in figure 3. The frequencies to the right of the mesh 
cut-off line are referred to as subgrid frequencies. If their production is 
permitted, they must alias back into the low-frequency range causing numerical 
error. This error can be severe enough to cause numerical instability. The 
standard way to cope with the subgrid scale generation is to include in the 
computing process some form of numerical dissipation which removes the sub- 
grid terms before any significant part of them cross the cut-off boundary. 
Notice that this is an arbitrary, numerical, error-control procedure that has 
nothing to do with any physical dissipation which occurs at much higher fre- 
quencies. 

The practical implementation of adding the numerical dissipation of the 
subgrid terms takes many forms. The process can be "hidden" in the differ- 
encing scheme. Such is the case for the various Lax-Wendroff types where the 
actual dissipative mechanism, which is provided by the fourth and higher 
even-order derivatives, is uncovered by inspecting the modified partial dif- 
ferential equation (e.g., Warming and Hyett, 197A, and Lerat, 1979). Upwind 
space-differencing schemes have the same property, which is again revealed by 
inspecting the modified partial-differential equation. Central differencing 
schemes for the first derivative of a space term are well known to be non- 
dissipative, so when these are used in shock-capturing algorithms, higher 
order dissipation terms are deliberately added to the computations (Von 
Neumann and Richtmyer, 1950, MacCormack and Baldwin, 1975, Warming and Beam, 
1976, Briley and McDonald, 1977, and Steger, 1978). From the arguments 
presented here, they are no better or worse than the forms which have no 
overt dissipation. All numerical schemes that capture shock waves with satis- 
factory accuracy have some numerical error, and its quantification is 
usually subjective and problem dependent. This situation can be attenuated 
to a certain extent by mesh clustering, but is usually worse for Navier- 
Stokes codes than it is for potential codes, simply because the meshes for 
Navier-Stokes computations are usually coarser. 

The above discussion presents one valid point of view for assessing 
shock-capturing techniques. However, it is not the only one. An alternative 
point of view is based on the theory of characteristics in supersonic flows. 
For example, the usual justification of upwind differencing in locally 
supersonic regions is not based on dissipation but on the fact that they can 
be made to approximate a local method of characteristics. The Lax-Wendroff 
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Figure 3.- Numerical dissipation of 
subgrid amplitudes . 
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methods also tend to approximate a local method of characteristics. In both 
cases, for all one-dimensional linear convective problems, a discontinuity 
can be solved exactly at a Courant number of one. Model problems, however, 
seldom occur in practical application. It is interesting to notice that 
central differenced first-derivative terms with deliberately added dissipation 
can be made to create a system that has the properties of upwind differencing. 

A second consequence of using a shock-capturing method is to create the 
problem of insuring the proper location and strength of the shock as it moves 
about in the mesh. Lax (1954 and 1973) has shown that this can be suitably 
approximated if the difference equations are locally conservative. The most 
common way of enforcing this condition is to cast the governing partial dif- 
ferential equations in conservation-law form, and then make sure the differ- 
ence scheme maintains this property. When such a technique is employed, a 
shock profile, represented, for example, by the pressure distribution, is 
"smeared’’ over a few mesh points, but, for many practical applications, the 
general position and strength are adequately represented. Many variations 
of shock-capturing methods exist which attempt to make the wave structure 
"crisper" and to eliminate overshooting of shock profiles. Our experience 
with numerical calculations which include boundary layer indicates that the 
details of shock smearing and overshoot are not of critical importance in 
determining the flow behavior along body surfaces. From this point of view, 
a wide variety of published methods are quite adequate for capturing shocks 
in Navier-Stokes codes. 


Numerical Techniques for Computing Turbulence Effects 

The problem of computing turbulence is much more difficult than that of 
capturing shocks. In fact, at the Reynolds numbers typical of transonic 
aerodynamic flows, no attempt is made to compute turbulence; rather, we try 
to approximate the effects of turbulence. The reason is, as in the case of 
shock capturing, the incapability of numerically resolving the full range of 
scale. However, in the case of turbulence, the problem is much more severe, 
since the scale to be resolved extends in all three space directions as well 
as in time. 

A plot of the longitudinal turbulence energy spectra for eight differ- 
ent types of flow is shown in figure 2. It is seen that energy-dissipating 
eddies (large k) are apparently independent of both Reynolds number and 
type of flow. Further, the form of the energy spectra in the inertial sub- 
range at high Reynolds number conforms to the Kolmogoroff spectrum law 
(k”?/!). This result is strictly experimental; no numerical simulation has 
yet produced real evidence of an inertial subrange in three dimensions. In 
order to accomplish this, one needs to provide a mesh that can support more 
than two orders of magnitude of frequency variation in all the three space 
dimensions. It is estimated that this will require a mesh with about (1024)^ 
grid points. For an incompressible flow that contains all of the modes, the 
calculation would need a total storage of about 7 x 10^ words using the most 
sophisticated numerical techniques. 
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The simple realities of computer resources force us to make one severe 
approximation and to accept one severe constraint in formulating our gov^ 
erning equations before we even start to consider the numerical methods. 

The approximation is the use of the Reynolds-averaged equations discussed 
earlier. This eliminates the need to resolve the small eddy motion, but 
introduces the problem of closure. The constraint is to permit extreme 
coordinate distortion in only one direction. This permits us to approximate 
viscous effects normal to very thin layers, but, at high Reynolds numbers, 
in that direction only. Probably the most important result of all this is 
that the computational processes that finally emerge have the capability 
of qualitatively simulating flows witii regions of separation and large-scale 
unsteady behavior. The crucial question, of course, is their reliability. 

We have yet to discuss the role numerics plays in computing turbulence 
effects. Two quite different issues are involved. One, the manner in which 
the subgrid scales are accounted for, and the other, the manner in which the 
turbulence model is implemented. The subgrid scales are constantly being 
generated by the large scale structure through the nonlinear wave interactions 
in the convective terms. The numerical control of the subgrid energy produc- 
tion is brought about by the addition of dissipation, either through the 
space derivative approximation or deliberately by additional terms. In either 
case the choice is arbitrary, except that it lie in the error band of the 
large scale resolution, and that it prevent the accumulation of energy in the 
highest frequencies supported by the mesh. The role of this form of dissipa- 
tion is often not clearly understood. It has absolutely nothing to do with 
physical viscosity at the scale that it is employed. Its detailed form is 
largely arbitrary, yet a solution would be physically incorrect if it were 
removed, since energy would then flow to subgrid levels and alias back into 
large-scale terms where it has no physical meaning. It is essential to the 
numerical simulation of the effects of turbulence, but it is not, in con- 
ventional terminology, part of the turbulence model, see figure 4. 

The second important role of numerics in Reynolds-averaged codes lies in 
the detailed coding of the turbulence model. The analytic forms of several 
models were given earlier in the section, Governing Equations. Unfortunately, 
these are not sufficient to describe the effect of a turbulence model on an 
actual calculation. The numerical effect of the complete model is the sum 
of all its parts, and this includes the grid clustering, the metric evaluations 
(see next section) , the internal logic controlling the local evaluation of 
parameters such as mixing length, and the choice of difference approximations. 
The "accuracy" of all this is difficult to evaluate since the conglomerate is 
the actual turbulence model and its fundamental basis is essentially empirical. 
The final judgement of the method is usually based on a comparison with some 
experiment, and the result may be good or bad depending on the choice of any 
one of the method constituents. 

Many variations of turbulence models have been tried on transonic flows 
with turbulent boundary layers. How well these compare with transonic wind- 
tunnel experiments is discussed in the next major section. 
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Figure 4 .- Modeling physics that cannot be computed. 
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Effect of Grid Choice on Numerical Stability 

The basic reason for choosing a nonuniform grid is to improve the 
accuracy of a numerical solution for a given number of mesh-points. There 
are two ways in which this is usually accomplished. One is to align, as 
closely as possible, a coordinate with a known or anticipated surface, such 
as a shock wave or body surface, in order to fit them more "naturally" into 
the mesh. The other is to cluster points in regions where there are rapid 
changes of gradients in order to reduce local truncation errors. As a 
corollary of the latter process, in order to conserve resources, points are 
often spread apart in regions where the curvature is small. Finite differ- 
ence, finite volume, and finite element methods all have these capabilities. 

The form chosen for a grid can have a profound effect on the solution 
process. By far the most important side effect of grid refinement on a 
numerical algorithm is its influence on numerical stability. As is very well 
known, the time step of explicit methods is mainly bounded by the size of the 
space interval, and this holds for nonequispaced as well as equispaced meshes. 
If a single time step is used for advancing the entire solution, an explicit 
method is generally limited by the smallest space interval in the mesh. This 
limitation can be seriously costly if the time step is forced to be very 
small compared to the time scales of motion that are of interest. In such 
cases, the algorithm is said to be stiff, and if the stiffness is caused by 
the fineness of a space interval in the grid, the algorithm is said to be 
mesh stiff. 


Codes using explicit numerical methods for the solution of the Reynolds- 
averaged Navier-Stokes equations can be extremely mesh stiff when they are 
used to study flows with thin boundary layers. This occurs when the grid is 
made to be very fine in the vicinity of the body in order to compute the 
viscous effects there. For example, a typical grid spacing normal to an air- 
foil surface can be in the order of 0.00001 chords for turbulent boundary- 
layer simulations at Reynolds numbers above 10&. Grid point clustering around 
shocks and leading and trailing edges can also be the cause of mesh stiffness. 

By far the most common way to avoid any form of stiffness is to use 
implicit, rather than explicit, algorithms. Almost all codes being used to 
analyze the compressible, Reynolds-averaged Navier-Stokes equations have some 
parts that represent an implicit numerical technique. The use of such tech- 
niques involves the solution of coupled sets of simultaneous equations. In 
finite difference codes, these simultaneous equations can usually be expressed 
as very sparse banded matrices that, in the great majority of cases, have 
tridiagonal structures. In fact, the numerical efficiency realizable from 
solving tridiagonal systems is so deeply embedded in finite difference methods 
for the problems we are discussing that it has greatly influenced, and at 
present even limits, the choice of grid topologies. This is discussed in the 
section. Effect of Grid Topologies on Computational Efficiency. 
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The Basic Difference Equations 

The following is a brief evaluation of the finite-difference techniaues 
currently being used to solve the types of problems in which we are interested 
norren^r consideration of the turbulence model. It does 

^ which compute the flow using different 

averaged J: llZlt efa 

werrj;v“oped';“““!d''r -<livldual cod'es 

nrohl^rr ?? used to solve specific problems. The particular choice of 

I rh ^ ^^k usually motivated by some experiment involving shock waves and 

te™?e ? "k“'’ «f separatfon. In ZJrIcfl 

terminology, these codes represent methods that range from fully explicit 

in re™^’ f ’ 1978)- The codes are usually written 

datl ™Lf "^hn^^ti ‘’T' “PPUP8 in series to prescribed 

nnpr^^ ’ Thus, there may be a convection operator followed by a diffusion 

r-o^eJar^r^L'^olf ^ one-dimensional 
X r so^uM^n f ^ one-dimensional y-operator to form the total 

to as ra^t^rL LLs techniques are also referred 

^ certain of the factors represent explicit 

methods and others implicit ones (MacCormack, 1978, and Shang, 1978). 

From a general point of view, at the present time, all of the codes 
sed to solve the Reynolds-averaged Navier-Stokes equations, and the methods 

cienc’^^'T^^"^-^^^"' potential for accuracy and effi- 

ciency of running time, although these can vary according to the^capabilities 

acce^Lblf^'' numerical methods they represent appeL to be 

hn.nS everywhere throughout the flow field except possibly at the 

matter which is again an individual responsibility. The codes 
g nerally at least first-order accurate in time. For high Reynolds 

to"*reach a steadv^st about 45 minutes of running time on a CDC 7600 

to reach a steady state, if one exists. This estimate is for codes that are 

of implicit. It varies, of course, depending upon thS number 

rhp^ H Reynolds number, and the angle of attack If 

the codes are fully explicit the running times can be much longer. 


Effect of Grid Topologies on Computational Efficiency 
expreI::d'J: the"L\r discussion our basic equation can be 



(27) 


where 

some 


A is a very large and 
combination of the flux 


very sparse nonlinear matrix that represents 
Jacobian, the grid construction, and the space 
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differencing. If the grid is chosen so that the physical space is mapped 
into a computational space that forms the inside of a rectangular box, and 
the boundary conditions are mapped onto the sides of the box, the matrix A 
becomes banded for most common choices of finite difference schemes. The 
typical form of A for second order finite-difference schemes is shown in 
equation (28) for a three-dimensional problem that is formulated in a computa- 
tional box. 



In this schematic structure, all matrix entries are zero except those repre- 
sented by the diagonal lines and each diagonal line represents a set of 5 x 5 
block matrices each of which is composed of a local flux Jacobian. Suppose 
the mesh coordinates are represented by x, y, and z and there are a total 
of Mx, My, and Mz points in each coordinate direction. In the particular 
case shown in equation (28), the data vector Q is so arranged that the x 
data is closely packed, nearby y data skip blocks of x, and nearby z 
data skip blocks of y. Of course, this arrangement is arbitrary and, by 
permuting the data base, the variables in any one direction can be closely 
packed at the expense of the other two. 
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One can view the steady state solution of equation (27) in two ways- 
one as the solution of the nonlinear system Q = A'lf, and tte other as the 
result of a converged time history of an unsteady process. The former 
requires the solution of a set of simultaneous equations having the form 
represented by A in equation (28)— which would have to be iterated because 
it is not linear. The latter would require the successive solution (with 
each time step) of a similar set of equations if the time-marching metho 
were fully implicit. 

Consider the prospect of carrying out either of these ' 

Although the matrix A is sparse and banded, notice that the half bandwidth 
is 5 X MX X My elements. A solution using simple Gaussian elimination would 
require about (5 x Mx x My) (5 x Mx x My x Mz) temporary storage locations to 
hold the information required to complete the backward sweep. Tins makes the 
solution of such a matrix by direct methods quite impractical on present day 
computers with even moderste mesh sizes. 

A common finite-difference technique used in the unsteady approach that 
overcomes the difficulty just discussed is to factor the t^e march process 
without changing the order of accuracy of the algorithm. There are s^ve^ 
ways for caryying this out, with differing accuracies and stabilities. They 
all have one thing in common which is to greatly reduce the temporary storage 
requirements for the implicit operation. Methods commonly referred to as 
factored fully-implicit lead to a set of three matrices representing block- 
tridiagonal equations that have to be solved in sequence. Each o t e 
matrices has the form shown in equation (29) . 



( 29 ) 
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Notice that this time the matrix is formed by large uncoupled diagonal blocks 
each one of which is tridiagonal in sub-blocks of 5 x 5 matrices. In such 
cases, each large diagonal block can be solved independently and requires a 
temporary storage of only 5 x 5 x Mp words, where p represents x, y, or z. 

The role of the topological-box computational space in all of this is to 
provide the banded structure of the matrices in equations (28) and (29). 

Zonal grids with inerfaces, overlapping meshes, and other forms of nonregular 
grid structures lead to A matrices that are not banded and tend to deviate 
from the tridiagonal structure. This can greatly increase the complexity of 
the computational algorithm or drive it to explicit (or even numerically 
unstable) forms. In either case, efficiency and code reliability can suffer. 
Many forms of the finite element approach will lead to the same difficulties 
for the same reason. The problem of generalizing mesh structures beyond com- 
putational boxes and keeping the codes that use them computationally reliable 
and efficient is one of the most pressing problems in finite difference 
developments in the 1980* s. 


A COMPARISON BETWEEN EXPERIMENTS AND CALCULATIONS 
OF TURBULENT TRANSONIC FLOWS 


The following material draws from the relatively young and limited body 
of computed results based on the Reynolds-averaged Navier-Stokes equations 
for transonic flows with strong viscous-inviscid interactions. We have taken 
this material from publications only from NASA Ames Research Center simply 
because most of the published work in this area has been carried out at this 
institution. 

First of all, consider some typical computed boundary-layer profiles for 
an attached flow. For example figure 5 shows a group of such profiles ahead 
of a shock wave on an 18% thick circular-arc airfoil. These are compared 
with the compressible form of the law of the wall. In the figure, u“‘“ repre- 
sents u normalized with the friction velocity (Deiwert, 1975). A simple 
mixing-length model, given by Launder and Spalding (1972), was used to 
describe the turbulent transport. All computed profiles have one grid point 
in the viscous sublayer. Notice that the log-law region is well represented 
by the grid-point distribution. This is generally the case of presently 
available Reynolds-averaged Navier-Stokes computations in attached boundary 
layers. Computed velocity values at x/c = 0.675 differ from the empirical 
log-law distribution because of flow separation just downstream of this 
location. 

For many practical uses, the turbulence modelling of attached turbulent 
boundary layers without shock-wave interaction is quite acceptable. The 
following discusses the status of the Navier-Stokes technology for turbulent, 
transonic simulations with emphasis on turbulence modelling for separated 
flows and on flow problems which are not feasible to solve with current 
simplified viscous-inviscid interaction approaches. This discussion deals 
with representative simulations in which flow fields may be steady, unsteady. 
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attached, or separated, both in two and three dimensions. Special problems 
such as "buffetting" flows, aileron buzz, and airfoil "stall" are considered. 
In addition, two types of turbulence simulations are presented, one in which 
turbulence models are used in a predictive mode and the other where these 
models are used in a postdictive mode. This presentation is motivated in 
order to stimulate systematic questioning of what research directions are 
needed for accelerating advances in better predictions of separated turbulent 
flows in aerodynamic applications. 


Axisymmetric Steady Flows 

A computation of normal shock-boundary layer interaction for an axisym- 
metric flow was carried out by Viegas and Horstman (1979). The tunnel geom- 
> experimental results, and several computations are shown in figure 6 . 

This represents an attempt to compare the merits of four different types of 
turbulence models at Reynolds number (5.5 x lo3) based on upstream boundary- 
layer thickness and low supersonic Mach number (1.44). The results for pres- 
sure distribution are essentially the same for all models. Differences are 
evident in calculation of skin friction, which depends, of course, on the 
slope of the boundary-layer profile at the surface. In fact, the particular 
algebraic model used showed a region of flow separation which did not appear 
in the other calculations. The above computers report that the most recent 
evaluation of the experiment indicates that flow does not separate. In light 
of the results shown in figure 7, a tentative conclusion can be drawn: 

This is probably representative of the accuracy one can expect from present 
forms of turbulence modelling and numerics. With regard to the algebraic 
model, the obvious question is: What details made the model used for Levy's 

results shown in figure 7 so superior to that used for the results in 
figure 6 ? 

As one looks into the details of more sensitive flow properties, one can 
anticipate further discrepancies. For example, the W-R model [equations 
(18) and (19) with 67 = 0.9] and the Jones-Launder (J-L) model (Jones and 
Launder, 1972) were used to compute the turbulent kinetic energy, v^/2. 

Measured and computed profiles of v2/2 are shown in figure 8 at various 
x-locations downstream of the shock wave located at x. The measured energy 
was determined from a measurement of U 2 and with the assumption that 

: U 2 : U 3 = 4 : 2 : 3. This assumption was observed to be reasonable for 
equilibrium boundary-layer flows at high subsonic Mach numbers (Acharaya, 1977). 

Computed Mach contours and the extent of separation region about an 
axisymmetric "bump" are shown in figure 9 , along with an infinite-fringe 
Interferogram and an oil-film visualization. A zero-equation model and the 
W-R model, respectively, predict shock locations 0.13 and 0.10 chord lengths 
downstream of the experimental location, which is at x/c ~ 0.66. Johnson 
and Horstman (1981) report that wall effects are negligible, and they believe 
the computational grid is sufficiently refined. Figure 9 also shows a 
surface oil-flow visualization indicating separation at x/c ~ 0 . 7 , and the 
experimental and computed locations of the uj = 0 line. 
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Figure 8.- Kinetic energy profiles 
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Figure 10.- Computed pressure distribution 





18% CIRCULAR-ARC AIRFOIL. Re = 11 x 10®. = 0.720. a = 0° 

(LEVY. 1978) 
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Figure 12.- Effect of tunnel-wall boundaries off-deslgn conditions. 


46 


18% CIRCULAR-ARC AIRFOIL, Re= 11 x 10®. = 0.783. a = 0° 

(LEVY, 1978) 


TUNNEL WALL 


STREAMLINE FROM FREE-FLIGHT SOLUTION 

1.00 1 


WALL COORDINATES y/c 


.96 

.92 


\ 


K- 


CONTOURED 
^ WALL REGION 


3 -2 -1 0 1 234567 
x/c 


PRESSURE DISTRIBUTION SKIN-FRICTION 



Figure 13.- Effect of tunnel-wall boundaries near-design conditions. 
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The wind-tunnel results shown in figures 10, 11, 12, 13, and 7 were 
determined from measurements made in a channel flow. Computations of such 
flows are known to be sensitive to inflow and outflow boundary conditions. 

In order to check this aspect of the problem, Coakley (private communication, 
1981) made some calculations in which the outflow pressure distribution was 
fixed at the experimental value. Some preliminary results are shown in 
figure 14. The results are encouraging. The shock wave is now oblique 
and the level of the trailing edge pressure plateau is matched quite closely. 
However, the shock position and skin friction are still parameter dependent 
and the investigation is continuing. 

We now turn to some results involving performance characteristics and 
shock-boundary layer interactions, but no shock- induced separation. 

The experimental and computed drag polars and lift curves for the GK I 
supercritical 11.5% thick airfoil (Garabedian and Korn, 1971) are shown in 
figure 15. The experimental data were taken with tunnel walls set a 6% and 
20,5% porosity (Kacprzynski , et al. , 1971). For 20.5% porosity, Melnik (1979) 
shows two sets of experimental data on the lift curve, uncorrected and cor- 
rected. According to him, the corrected data represent free-flight conditions 
(see also Morky and Ohman, 1980). There are two sets of computed results. 
Deiwert (private communication, 1977) has solved the Reynolds-averaged Navier- 
Stokes equations with free-flight boundary conditions and an algebraic model 
without relaxation. Melnik (1979) has used the ’’full" viscous-inviscid inter- 
action theory. He has matched the lift coefficient with the experiments and 
applied a small Mach number shift of M = -0.005 to obtain agreement with 
the experimental shock position. The lift curve shows that both the viscous 
effects and the wind-tunnel interference effects are important. Drag values 
of both computations differ from the measured values. These computations 
again indicate that proper boundary conditions are required for taking into 
account wind-tunnel wall-interference effects. 


Two-Dimensional Unsteady Flows 


An interesting set of experiments (McDevitt, 1976) and calculations 
(Levy, 1978) have been carried out for an 18%-thick biconvex airfoil at zero- 
degree incidence. Both experiments and calculations showed a region of 


”buffetting" or self-excited, oscillating flow in the Mach number range 
between 0.72 and 0.79 for a Reynolds number around 11 x 10^. 


The experiment was conducted using a wind tunnel in which the upper and 
lower walls were contoured as mentioned in the preceding subsection. The 
calculations used slip-flow boundary conditions along surfaces that matched 
these contours. The effect of turbulence was approximated by an algebraic 
eddy viscosity model similar to that used by Deiwert (1977). This zonal 
model changed form in various regions bounded by the separation location, the 
location of reattachment of the separated streamline to the surface stream- 
line and the edge of the boundary layer. Unfortunately, the sensitivity of 
the solution to the model is an unknown. 
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Figure 15,- Lift coefficients and drag polar for GK I airfoil- 
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Figures 16 to 19 show a comparison between the experiment and computed 
results. Figure 16 identifies the experimental Reynolds number and Mach 
number domains within which there are three distinctively different types of 
flow. The three types were reproduced by computations made at Mach numbers 
of 0.720, 0.754, and 0.783, and a chord Reynolds number of 11 x lOo. At 
— 0.720, the flow is steady and flow separation occurs near the trailing 
edge of the airfoil. At M«> = 0.754, there is unsteady periodic oscillation 
in shock— wave location and intensity; and the flow alternates between trailing- 
edge and shock- induced separation and is quite different on the upper and 
lower surface at any given time. At Mco = 0.783, a shock wave induces 
boundary-layer separation at its base and the flow is relatively steady, 
except in the separated region. 

Surface pressure comparison is demonstrated between computations and 
experiments for the above three different conditions in figure 17. The 
vertical bars on the experimental data represent maximum and minimum values 
of fluctuations about mean. The range of computed fluctuations about the 
mean computed values is denoted by the shaded area. The steady flow regions 
at Moo = 0.720 and 0.783 have been discussed in the previous subsection. 

The unsteady flow at Moo = 0.754 is qualitatively very well predicted, 
but quantitative comparison is poor, except for the mean values of pressure 
over the forward half of the airfoil (figure 17). This is further supported 
by figure 18 which shows surface-pressure time histories. Here, the instan- 
taneous pressure oscillations are given about the mean pressure, normalized 
by the wind-tunnel total pressure. The computed and measured, reduced fre- 
quency of these oscillations are, respectively, 0.40 and 0.49. However, the 
amplitude of oscillations is quite different. For this case, the shock-wave 
shapes from shadowgraphs are compared with computed Mach number contours in 
figure 19 where the phase has been arbitrarily adjusted (Marvin et al., 1980). 
For another problem, namely, a 14%— thick biconvex airfoil at Re — 7 x 10 
and Moo = 0.83, the computed unsteady lift forces and pitching moments are 
compared with those for Moo = 0.85 in figure 20 (Levy, pr i v a t a rnimmini nr i on ., 

1981) . 

It is not at all surprising that Reynolds-averaged Navier-Stokes equa- 
tions are capable of simulating unsteady flows when the computational time- 
step is small compared to the period of resolvable flow motion which is of 
interest, but much larger than the high-frequency, small-scale fluctuations 
which have been averaged out of these equations (see earlier section. 

Governing Equations). The question of how high the resolvable frequency 
could be relative to the mean frequency of turbulence eddies is addressed by 
Chapman (1979). 

Another unsteady phenomenon, this time associated with a moving boundary, 
is represented by the performance characteristics of the aileron of a P-80 
(l.e., F-80) aircraft. This flow has been simulated by Steger and Bailey 
(1980) using the algebraic eddy viscosity model and the second-order thin- 
shear— layer approximation described in the section. Governing Equations. The 
turbulence model w..is applied from the leading— edge of the airfoil. The P-80 
airfoil section is an NACA 652^-213 with a = 0.5. The aileron buzz is a 



STEADY FLOW, UNSTEADY FLOW, STEADY FLOW. 
TRAILING EDGE OSCILLATORY SHOCK-INDUCED 

SEPARATION SEPARATION SEPARATION 



Figure 16.- Experimental flow domains for 18% circular -arc airfoil. 
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18% CIRCULAR-ARC AIRFOIL, Re = 11 x 10 





18% CIRCULAR- ARC AIRFOIL, Re = 11x 10®, M^=0.76,a=0° 
(SEEGMILLERet al„ 1978) 


UPPER SURFACE 

x/c =0.50 x/c = 0.775 

.2 n .2 r 



LOWER SURFACE 



0 9 18 27 0 9 18 27 

CHORDS TRAVELED CHORDS TRAVELED 


Figure 18.- "Buffeting" flow, surface-pressure time histories. 
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one-degree-of“freedom flutter problem (Erlkson and Stephenson, 1947). The 
interrupted, inviscid shock-wave motion (e.g., Tijdeinan, 1980) causes a phase 
shift in the response of the hinge moment to the aileron movement. In the 
experiment, at Moo = 0.82 and a = “1°, the aileron, when freed at an angle 
near zero, would buzz with amplitude and frequency as indicated in figure 21. 
In the simulation, it would not buzz under these conditions. But if it was 
initially deflected to 4°, it would, on being released, buzz as shown in the 
figure. The computed and measured frequency are, respectively, 22.2 Hz and 
21.2 Hz. Further, the computed and measured deflection of the aileron are, 
respectively, -1,1 ± 11.1 and -3 ± 9.2 deg. Similar calculations are made at 
different airfoil angles of attack to predict the measured buzz boundary. 

Figure 22 shows results for an unsteady transonic flow over an NACA 
64A010 airfoil, which is oscillating about its one-quarter chord with a 
reduced frequency of 0.2, based on one-half chord. Chyu et al. (1981) 
obtained these results with the same CDC 7600 computer code used for the buzz 
study discussed above. The computations were done in a coordinate system 
fixed to and moving with the airfoil, but stationary at the open boundaries. 
This involved generation of a grid system for each time step. The above 
investigators report no flow separation. Computed and measured surface pres- 
sure distributions are shown only for one-half cycle of an oscillation, as the 
airfoil angle varies from 1 deg to -1 deg. Notice that the computed and mea- 
sured results agree much better downstream from the shock wave than upstream 
of the shock. Figure 23 shows computed and measured shock-wave locus on the 
upper surface of the airfoil. 

Recently, '*stall" boundary of the GK I airfoil has been predicted by Levy 
and Bailey (1981). The Illiac IV computer code was the same as that used on 
the buzz study just discussed. Figure 24 shows computed and measured unsteady 
flow boundaries^ and computed Mach contours. This figure shows much better 
agreement between experiment and calculations at the high-Mach-number , low- 
lift range than they do on the low-Mach-number , high-lift side. The latter 
represents a case where a turbulence model has been pushed far beyond its 
limits. Notice the Mach contour plots in figure 24 at two different free- 
stream Mach numbers. In the low-Mach-number case, there is shock- induced, 
turbulent separation bubble. Whether in an experiment there is a transitional 
bubble ahead of the shock wave or below it, remains to be determined. In the 
high-Mach number case, there is again shock- induced separation which extends 
beyond the trailing-edge of the airfoil. 


Three-Dimensional Steady Flows 

In their present forms, most Reynolds-averaged Navier-Stokes codes for 
two-dimensional flows take rather lengthy, 0.75 to 3.5 hours on a CDC 7600, 
run times for grids of 4 to 10 thousand points to reach a steady state or the 
onset of a periodic flow. Three-dimensional flow simulations on such computers 
are, therefore, not common. On the so-called class VI computers, such as the 
ILLIAC IV, however, some three-dimensional studies with moderate resolution are 
practical at a research level. We conclude with a brief discussion of two of 
these investigations. 


^Investigators of these boundaries have called them buffet boundaries, although 
there was no aeroelastic response of the airfoil to aerodynamic excitation 
arising from unsteady separated flow (Fung, 1955). 
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COMPUTATION OF AILERON BUZZ BOUNDARY 



• BUZZ 



ANGLE OF ATTACK, deg 



Figure 21.- A performance characteristic of the control surface 

of P80 aircraft. 
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OSCILLATING NACA 64A010 AIRFOIL, Re= 1.2x 10^, = 0.8. k - 0.2 

(CHYUet al., 1981) 



° lS experiment 
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OSCILLATING NACA 64A010 AIRFOIL, Re = 1.2 x 10^, = 0.8, k = 0.2 

(CHYUet al., 1981) 
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INVISCID I COMPUTATION 
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Figure 23.- Shock wave locus on upper surface. 
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UNSTEADY FLOW BOUNDARIES 



Figure 24.- A performance characteristic in maneuver of Gar abed ian -Korn 

airfoil at Re = 21 x 10^. 




Surface pressure isobars for a subcritical unseparated flow over a 45° 
swept, 10%-thick circular-arc airfoil at a zero incidence and spanning a 
tunnel are shown in figure 25. The second-order thin-shear-layer approxima- 
tion was used with the two-equation W-R turbulence model, and only the upper 
half of the flow field was computed. 

Notice that the computational and experimental Mach numbers are slightly 
different. Bertelrud et al. (1980) have explained the difference between the 
Mach numbers as follows: The reference Mach number and pressure values for 

the experimental isobars were obtained at a location nearly one chord length 
ahead of the wing leading edge at the left wall of the channel. The computa- 
tional boundary was located at 3.5 chord lengths ahead of same leading edge 
of the wing. Therefore, the computed state at the measuring location did not 
correspond to the measured state at that location. Figure 25 shows a compari- 
son of measured and computed pressure distributions at three spanwise loca- 
tions on the wing surface, and it gives the Mach number sensitivity. 

Simulations of three-dimensional boat tail afterbody flow fields have 
been obtained by Deiwert (1980) with the second-order thin-shear-layer approxi- 
mation and the same algebraic turbulence model used in the buzz study dis- 
cussed above. In figure 26, surface pressure distributions are shown for a 
boattail model used by Shrewsbury (1968). The experimental data are shown in 
the insert by the triangles, squares and circles corresponding to windward, 
lateral, and leeward positions. The corresponding computed results are shown 
by dashed, dotted, and solid lines. The junction of the forebody and after- 
body of the above boattail model is sharp. Deiwert (1980) has reported some 
sensitivity of the computed results to the grid spacing in the vicinity of this 
junction (figure 27). Figure 28 shows computed results. The upper part is sur- 
face pressure topology and the lower one is a limiting surface flow pattern 
(surface shear directions) which approximates a surface oil-flow pattern. 

The symbols S and R, respectively, stand for flow separation and flow 
reattachment; and the subscripts S and N, respectively, denote a saddle-point 
and a node-point. Downstream of the circumferential line Sg Sjq, the flow is 
separated. Downstream of the circumferential line Rs the flow is 

attached. The direction flow is from Sg to and from R^ to Rg. Such 
details are available from present Navier-Stokes technology, and they are of 
considerable use towards a better understanding of complex flow fields and 
towards providing internal consistency checks for simulations. 


CONCLUDING REMARKS 


The Navier-Stokes technology is currently under vigorous development. 

It has opened new possibilities of simulating unsteady, separated, turbulent, 
compressible flows that were not accessible five years ago. In this paper we 
have presented the state of the art, as we envision it. The primary utility 
of this technology is in applications where the present viscous— inviscid inter- 
action computations fail, and this generally occurs in simulating separated 
flows that are nominally two-dimensional and unsteady, or three-dimensional 
steady or unsteady. There is little doubt that the Navier-Stokes technology 
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Re = 2.9 X 10®, = 0.9, a = 6° 

(DEIWERT, 1980) 

ODA EXPERIMENT 

(Shrewsbury: 

1968) 

— — JCOMPUTATION 



.2 .6 1.0 1.4 


X + 0.638 

Figure 26.- Comparison of surface-pressure distribution 
on boattail afterbody. 
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= 0-9. Rod = 2.9 X 10®. a = 
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AFTERBODY GRIDS 




AFTERBODY PRESSURE DISTRIBUTION 


O EXPERIMENT (SHREWSBURY, 1968) 

COMPUTATION, UNIFORM GRID 

COMPUTATION, CLUSTERED GRID 

(DEIWERT, 1980) 


Figure 27.- Influence of forebody/afterbody juncture. 
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developers; the oeestlon Is ho„ 


i° substantial advances have been made In computer 

speed and mmory. Further, our ability to compute flows in rathe^ cZier 
? d® greatly Improved. It is becoming increasingly clear that 

turbulence modelling in the regions of separation, which are not "small " is 
the weakest part of the Navier-Stokes technology. In facr it L 
becoming the primary pacing item for Reynolds-averaged Navier-Stokes. ^ 

prove the capability of 

simulating turbulent flows with mild separation. Many turbuLnce mLels 
were made to work on Isolated experiments. In fact this approach was used to 
models^ ^ ^ numerical techniques and a few empirical cLstants in the 

models. However, very little has been done to establish the reliabilitfof 

must remember that the same turbulence "model" can give diffLent 

methin r different codes with the same or different n^er^cff 

methods. This is due principally to lack of grid-refinement stud^^s! 

Numerical simulations of the Reynolds-averaged Navier-Stokes equations 
are, in general, predictive for attached boundary layers. Zero-equatio^ 

^uff he : engineering analysis'of these "lows! Jut they 

and fin with caution when used to approximate separated flows 

and nows with strong curvature effects. Simplicity of zero-equation models 
mnre'^e^ adjustment for separated flows; complex models, which contain 

this tiL^ constants, need less adjustment. From the results available at 
tf Low ^ero-equation models are judged, there is no clear evident 

to show that one- or two-equation, first-order models are much better. 

Is d..e"! problems in constructing models for external separated flows 

ilt f, ?? the behavior of turbuleLHn 

instaLr^tW^?^^^'^’ Johnston, 1980). We do know, for 

lenL L’ separated flows normal stresses are anisotropic and turbu- 

ence structure xs not in equilibrium. Relaxation procedures and transport 
equations for turbulent scales can take into consideration some of the Mstory 
effects, namely, the nonequilibrium nature of turbulence; but the Reynolds ^ 

f? TLr IS modelled to respond instantly to changes in mean strain 

field [equation (10)]. The first-order (eddy-viscosity) models, therefore 
can be truly predictive only for flows in which turbulence is nLrly in !Lal- 
equl xbrium 6r for self-preserving flows. The second-order (s!!!ssLquai!on) 
models are required for nonequilibrium flows. This is illusLated b!!o! 

Consider a distortion of a flow field of fully developed, homogeneous 

of plane strain (figure 29). This experiment acts 
are aTselt ^hf separated flows when near-surface effects 

LL w ^ conditioned through screens, and it becomes parallel 

turbulence to become anisotropic. A measure of anisotropy is plLted as the 
ordinate, the lower portion of figure 29. At some distaLe doLstrel? tL 
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two-equation MODEL 

1 [ 1 1 1 — 1 1 

CONST RATE PARALLEL 

OF STRAIN "" FLOW 

.8 - 



X, m 


Figure 29.- Normally strained homogeneous flow (Wilcox and Rubesln, 1980) 




strain is removed and the fluid returns to parallel flow. The measurements 
of Tucker and Reynolds (1968) are compared with computed results of Wilcox 
and Rubesin (1980). The computation with a second-order (Reynolds stress) 
model gives a better agreement with the measured values than that of the 
^^J^®t-order (W-R) model. Although Wilcox and Rubesin modified equation (10) 
to remove the alignment of the Reynolds stress tensor and the mean rates of 

model, the predicted return to Isotropy is abrupt when the 
strain is removed. This kind of behavior is brought about by the shortcoming 
of the first-order models as explained above. 

The above example illustrates two points. First, eddy— viscosity models 
can probably never be completely predictive for separated flows. Some details 
of the structure will most certainly be lacking. More sophisticated models 
wij.1 pick up more of the details, but for stringent requirements they, too, 
may fail. The second point, and by far the most important one, is that it 
is probably possible to predict the gross behavior of a flow even when certain 
of the details are not well represented or even missing altogether. The most 
meaningful test of whether or not these calculations have useful information 
is whether or not they are used. 

There are two schools of thought about modelling turbulence (Lumley, 
1978). Some believe that under certain circumstances, rational second-order 
(or invariant) modelling can be developed for general computation procedures. 
They consider this approach may at least provide a guide for the construction 
of the more empirical models. Others believe the structure of turbulence to 
be so complex that a search for universal closures is probably in vain. 

They believe that practical computations will require empirical techniques 
developed for particular flow topology. As for the current efforts in com- 
puting turbulence flows for industrial needs, Liepmann (1979) has presented 
an adversely critical opinion. 

There are probably five different parallel avenues of turbulent separated 
flow research: (1) Different turbulence models are applied to the same 

geometrical flow problem in order to determine which one is the best; 

^he same model, without any change in its form or in its empirical con- 
stants, is applied to different geometrical flow problems so that its breadth 
of application can be determined; (3) for a given form of a model, a computer 
optimization is carried out to obtain the best set of model parameters rela- 
tive to an available set of experiments; (4) for a specific flow problem, a 
determination of the range of flow parameters is carried out for which a 
given model with its empirical constants is valid; (5) a model is developed 
for a particular flow problem based on a detailed experimental characteriza- 
tion of this flow. As demonstrated earlier, elements of avenues (1) and (2) 
are already being pursued; remaining avenues need to be pursued. 

Computational aerodynamics is probably going to depend more on experi- 
mental inputs and checks and less on the solutions of the Navier-Stokes 
equations for developing turbulence models. Therefore, experimentalists 
should be requested to document well the experiments they conduct during their 
quest for understanding turbulence in separated flows. Through that under- 
standing, better turbulence models may result at least for these flows, and 


69 



this should increase the utility of the Navier-Stokes technology. Both 
experimental and theoretical investigators need to work together to advance 
the state of the art of turbulence models for separated flows. It is hoped 
that efforts will be devoted to extensive testing of these models on a variety 
of experiments without modifications to the basic coding. In addition, 
repeated grid-refinement studies are required to demonstrate that a turbulent 
numerical simulation tends to be independent of numerics. 

In the 1980’ s, the complex three-dimensional geometries will require 
component-adaptive or zonal methods. These procedures, along with limited 
availability of computer speed and memory, will guide the Navier-Stokes 
technology towards a viscous-inviscid interaction approach, which probably 
will consist of matching the Reynolds-averaged Navier-Stokes solutions next 
to a body surface with either Euler or potential flow solutions away from 
the surface. 

If the above efforts prove to work then not only capability but reli- 
ability would be established. At this point the Navier-Stokes technology 
will come of age. 

In summary, the state of the art of viscous transonic aerodynamics is 
presented in a Venn diagram shora in figure 30. At present, transonic, 
attached, two-dimensional, steady and fully turbulent flows can be routinely 
predicted. Extensive efforts are being made to predict both steady and 
unsteady, two-dimensional fully turbulent separated flows. Already promising 
starts have been made to simulate steady three-dimensional flows, either 
attached or separated. However, much remains to be done for laminar- 
transitional-turbulent flows. Further, there is negligible progress in meet- 
ing the final objective of predicting unsteady, three-dimensional, separated, 
and laminar- transitional- turbulent flows. 
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Figure 30. - Viscous transonic aerodynamic 
the state of the art. 
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